subroutine initial
    use commondata
    implicit none
    real*8 :: rho, u, v, p, c, b    ! b--vortex intensity
    real*8 :: xbar, ybar, r2    !  relative position in isentropic vortex problem
    real*8 :: btemp
    
    open(101, file="ini.txt")
    read(101,*) CFL             !  CFL Number
    read(101,*) FLG_prob        !  Problem:    1__Rayleigh Taylor;    2__Double Mach Reflection;    3__Forward Step;
    read(101,*)                 !  Problem:    4__Isentropic Vortex;  11--17__2D Riemann_1--7
    read(101,*) FLG_split       !  Splitting:    !!1__local-L-F;    !!2__global-L-F;    3__Steger-Warming;    4__VanLeer
    read(101,*) FLG_order       !  Order:   5__5th order;   7__7th order
    read(101,*) FLG_recon       !  Reconstrution:    1__WENO;    2__WENO-z;    3__WENO-zn
    close(101)

    
    if(FLG_prob==1) then    ! Rayleigh Taylor
        im = 120
        jn = 480
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 5.d0/3.d0
        xlong = 0.25d0;    ylong = 1.d0;    t_end = 1.95d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "Rayleigh Taylor"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( ys(i,j) <= 0.5d0 ) then
                    rho = 2.d0;    p = 2.d0*ys(i,j) + 1.d0;    c = sqrt(gamma*p/rho)
                    u = 0.d0;      v = -0.025d0*c*cos( 8.d0*pi*xr(i,j) )
                else
                    rho = 1.d0;    p = ys(i,j) + 3.d0/2.d0;    c = sqrt(gamma*p/rho)
                    u = 0.d0;      v = -0.025d0*c*cos( 8.d0*pi*xr(i,j) )
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==2) then    ! Double Mach Reflection
        im = 960
        jn = 240
        allocate( xr(-4:im+4,-4:jn+4), ys(-4:im+4,-4:jn+4) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 4.d0;    ylong = 1.d0;    t_end = 0.2d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "Double Mach Reflection"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j) <= (1.d0/6.d0 + ys(i,j)/sqrt(3.d0)) ) then
                    rho = 8.d0
                    u = 8.25d0*sqrt(3.d0)*0.5d0
                    v = -8.25d0*0.5d0
                    p = 116.5d0
                    c = sqrt(gamma*p/rho)
                else
                    rho = 1.4d0
                    u = 0.d0
                    v = 0.d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==3) then    ! Forward Step
        im = 600
        jn = 200
        ifs = im/5
        jfs = jn/5
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 3.d0;    ylong = 1.d0;    t_end = 4.d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "Forward Step"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        rho = 1.4d0
        u = 3.d0
        v = 0.d0
        p = 1.d0
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==4) then    ! Isentropic Vortex
        im = 64
        jn = 64
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0;    b = 5.d0
        xlong = 10.d0;    ylong = 10.d0;    t_end = 100.d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "Forward Step"'
        write(100,*) 'variables=', '"x","y","rho","u","v","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                xbar = xr(i,j) - 5.d0
                ybar = ys(i,j) - 5.d0
                r2 = xbar**2 + ybar**2
                btemp = -(gamma-1.d0)*b**2*exp(1.d0-r2)/(8.d0*gamma*pi**2)
                rho = (1.d0+btemp)**(1.d0/(gamma-1.d0))
                u = 1.d0 - b/(2.d0*pi)*exp(0.5d0*(1.d0-r2))*ybar
                v = 1.d0 + b/(2.d0*pi)*exp(0.5d0*(1.d0-r2))*xbar
                p = rho**gamma
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, u, v, p
            end do
        end do
        
    else if(FLG_prob==11) then    ! 2D Riemann  --1
        im = 400
        jn = 400
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.8d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -1"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)<0.8d0 .and. ys(i,j)<0.8d0 ) then
                    rho = 0.138d0
                    u = 1.206d0
                    v = 1.206d0
                    p = 0.029d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.8d0 .and. ys(i,j)>=0.8d0 ) then
                    rho = 1.5d0
                    u = 0.d0
                    v = 0.d0
                    p = 1.5d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<0.8d0 .and. ys(i,j)>=0.8d0 ) then
                    rho = 0.5323d0
                    u = 1.206d0
                    v = 0.d0
                    p = 0.3d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.8d0 .and. ys(i,j)<0.8d0 ) then
                    rho = 0.5323d0
                    u = 0.d0
                    v = 1.206d0
                    p = 0.3d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==12) then    ! 2D Riemann  --2
        im = 2048
        jn = 2048
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.25d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -2"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)<0.5d0 .and. ys(i,j)<0.5d0 ) then
                    rho = 0.8d0
                    u = 0.d0
                    v = 0.d0
                    p = 1.0d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 0.5313d0
                    u = 0.d0
                    v = 0.d0
                    p = 0.4d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 1.d0
                    u = 0.7276d0
                    v = 0.d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)<0.5d0 ) then
                    rho = 1.d0
                    u = 0.d0
                    v = 0.7276d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==13) then    ! 2D Riemann  --3
        im = 400
        jn = 400
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.23d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -3"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)<0.5d0 .and. ys(i,j)<0.5d0 ) then
                    rho = 1.d0
                    u = 0.75d0
                    v = 0.5d0
                    p = 1.0d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 1.d0
                    u = -0.75d0
                    v = -0.5d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 2.d0
                    u = -0.75d0
                    v = 0.5d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)<0.5d0 ) then
                    rho = 3.d0
                    u = 0.75d0
                    v = -0.5d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==14) then    ! 2D Riemann  --4
        im = 400
        jn = 400
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.2d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -3"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)>=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 1.d0
                    u = 0.d0
                    v = 0.d0
                    p = 1.0d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 0.5197d0
                    u = -0.7259d0
                    v = 0.d0
                    p = 0.4d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 1.d0
                    u = -0.7259d0
                    v = -0.7259d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 0.5197d0
                    u = 0.d0
                    v = -0.7259d0
                    p = 0.4d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==15) then    ! 2D Riemann  --5
        im = 400
        jn = 400
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.25d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -3"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)>=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 1.5d0
                    u = -0.75d0
                    v = -0.5d0
                    p = 1.5d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 2.d0
                    u = -0.75d0
                    v = 0.5d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 1.d0
                    u = 0.75d0
                    v = 0.5d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 3.d0
                    u = 0.75d0
                    v = 0.5d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==16) then    ! 2D Riemann  --6
        im = 400
        jn = 400
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.25d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -3"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)>=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 0.5197d0
                    u = 0.1d0
                    v = 0.1d0
                    p = 0.4d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 1.d0
                    u = 0.6259d0
                    v = 0.1d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 0.8d0
                    u = 0.1d0
                    v = 0.1d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 1.d0
                    u = 0.d0
                    v = 1.206d0
                    p = 0.3d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    else if(FLG_prob==17) then    ! 2D Riemann  --7
        im = 400
        jn = 400
        allocate( xr(0:im,0:jn), ys(0:im,0:jn) )
        allocate( var(kn+1,-4:im+4,-4:jn+4), v0(kn,-4:im+4,-4:jn+4), v1(kn,-4:im+4,-4:jn+4), vtemp(kn,-4:im+4,-4:jn+4) )
        allocate( fxp(kn,-4:im+4,-4:jn+4), fxm(kn,-4:im+4,-4:jn+4), fyp(kn,-4:im+4,-4:jn+4), fym(kn,-4:im+4,-4:jn+4) )
        allocate( fx(kn,-1:im,-1:jn), fy(kn,-1:im,-1:jn) )
        
        gamma = 1.4d0
        xlong = 1.d0;    ylong = 1.d0;    t_end = 0.25d0
        dx = xlong/real(im);    dy = ylong/real(jn)
        open(100, file='initial.dat')
        write(100,*) 'title = "2D Riemann -3"'
        write(100,*) 'variables=', '"x","y","rho","p"'
        write(100,*) 'zone t= "', tn, '",i=',jn+1, ',j=',im+1, ',f=point'
        do i = 0, im
            do j = 0, jn
                xr(i,j) = real(i)*dx
                ys(i,j) = real(j)*dy
                if( xr(i,j)>=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 0.5313d0
                    u = 0.1d0
                    v = 0.1d0
                    p = 0.1d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)>=0.5d0 ) then
                    rho = 1.0222d0
                    u = -0.6179d0
                    v = 0.d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)<=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 1.d0
                    u = 0.1d0
                    v = 0.1d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                else if( xr(i,j)>=0.5d0 .and. ys(i,j)<=0.5d0 ) then
                    rho = 1.d0
                    u = 0.1d0
                    v = 0.8276d0
                    p = 1.d0
                    c = sqrt(gamma*p/rho)
                end if
                v1(1,i,j) = rho
                v1(2,i,j) = rho*u
                v1(3,i,j) = rho*v
                v1(4,i,j) = p/(gamma-1.d0) + 0.5d0*rho*( u**2 + v**2 )
                write(100,*) xr(i,j), ys(i,j), rho, p
            end do
        end do
        
    end if
    
    close(100)
end subroutine initial  